Astronomy <fc Astrophysics manuscript no. physical'params'jets 


February 2, 2008 


(DOI: will be inserted by hand later) 





Which physical parameters can be inferred 
from the emission variability of relativistic jets? 

P. Mimica 1 , M.A. Aloy 1 ' 2 , E.Miiller 1 and W. Brinkmann 3 

1 Max-Planck-Institut fur Astrophysik, Postfach 1312, D-85741 Garching, Germany 

2 Departamento de Astronomi'a y Astrofi'sica, Universidad de Valencia, 46100 Burjassot, Spain 

3 Max-Planck-Institut fur extraterrestrische Physik, Postfach 1603, D-85740 Garching, Germany 

Received ?; accepted ? 

Abstract. We present results of a detailed numerical study and theoretical analysis of the dynamics of internal 
shocks in relativistic jets and the non-thermal flares associated with these shocks. In our model internal shocks 
result from collisions of density inhomogeneities (shells) in relativistic jet flows. We find that the merged shell 
resulting from the inelastic collision of shells has a complicated internal structure due to the non-linear dynamics 
of the interaction. Furthermore, the instantaneous efficiency for converting kinetic energy into thermal energy is 
found to be almost twice as high as theoretically expected during the period of significant emission. The Lorentz 
factors of the internal shocks are correlated with the initial inertial masses of the shells. 

Because of the complexity of the non-linear evolution the merged shell becomes very inhomogeneous and simple 
one-zone models are inadequate to extract physical parameters of the emitting region from the resulting light 
curves. In order to improve on these one-zone approximations, we propose a novel way of analyzing the space-time 
properties of the emission. Based on these properties we construct an analytic model of non-thermal flares which 
can be used to constrain some (unobservable) physical parameters of the internal shocks. These are the ratio of the 
Lorentz factors between the forward and the reverse shock (caused by the shell collision), and the shell crossing 
times of these shocks. The analytic model is validated by applying it to the synthetic light curves computed from 
our models. It can equally well be applied to observations. 

Key words, galaxies: BL Lac objects: general; X-rays: general — radiation mechanisms: non-thermal; acceleration 
of particles; method:numerical:hydrodynamics 



1. Introduction 

BL Lac objects are a class of active galactic nuclei 
(AGN) which show the most rapid variability of all AGNs. 
In the X-ray frequency range they display flares whose 
duration in the observer frame is usually of the or- 
der of one day ifMaraschi et al. I 1999, Takaha shi et al. I 
2000, Kata oEa et al. I 2001). as further confirmed by the 
longterm observation of Mrk 501, PKS 2155-304, and 
Mrk421 IjTanihata et al. l 2001. lTanihatal 200lL With the 
highly improved sensitivity of XMM — Newton, variabil- 
ity on time scales down to a few kiloseconds could be stud- 
ied in Mrk 421 ( Brinkm ann et al. 2001, 2003) and only re- 
cently, from an XMM — iVeiufonobservation of Mrk 421 
in a very high state, the spectral evolution of the ob- 
ject down to time scales of > 100 sec could be followed 
IjBrinkmann et al. 1 2005) . 

The internal shock scenario IjRees fe Meszarosl 1994) 
is usually invoked in order to explain the variability 
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of blazars flSpada et al~| 2001, |Bicknell fc Wagner| 2002). 
These shocks are typically assumed to arise due to the 
collisions of density inhomogeneities within the relativis- 
tic blazar jets. Recently, one-dimensional l|Kino et al. I 
2004) and two-dimensional (IMimica et al. I 2004; there- 
after MAMB04) simulations of such internal shocks have 
been performed both showing that the evolution of inter- 
nal shocks is considerably more complicated than what has 
been concluded from approximate analytic considerations. 
As a result of the complexity of the non-linear evolution 
arising after a two-shell interaction, the merged structure 
becomes very inhomogeneous and simple one-zone models 
are inadequate to extract physical parameters of the emit- 
ting region from the resulting light curves. In particular, 
we will show that the rest-mass density of the emitting 
shell can hardly be inferred from fits to the synthetic light 
curves. Assuming an approximately uniform proportion- 
ality between the rest-mass energy and the magnetic field 
energy the magnetic field strength might hardly be in- 
ferred from a flare light curve, too. Should the evolution 
of real shell interactions be as complex, the parameters 
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obtained from one-zone model fits to the spectra will be 
equally inaccurate as those obtained from one-zone model 
fits of the light curves. 

In this paper we both extend our simulations of inter- 
nal shocks and perform a systematic study of the influence 
of the fluid properties on the formation of flares, on their 
duration and temporal profiles, and develop an analytic 
model linking the fluid properties to the observed flare 
properties. One of the advantages of our analytic model 
is that it can be validated by comparing its predictions to 
the actual physical conditions obtained from our relativis- 
tic hydrodynamic simulations. This validation will probe 
to be very important considering the disparity of physical 
values predicted by other analytic models for single flares 
in certain blazars (see §0). 

In section §[2 we discuss the internal shock scenario as 
a model which explains the flare emission and we summa- 
rize our previous findings (MAMB04) regarding the hy- 
drodynamic evolution of internal shocks. The numerical 
simulations are described in Section §[21 Finally, in §0]we 
present an analytic model which can be used to infer phys- 
ical parameters of the fluid flow from either synthetic or 
observed flare light curves (§E] and[BJ . 

2. Internal shock scenario 

The internal shock scenario assumes the existence of blobs 
of matter moving with different velocities along the jet, 
presumably emitted by the intermittently working central 
engine. Two blobs will interact after a certain time de- 
pending on their initial relative velocity. Once the interac- 
tion of the shells starts two internal shocks form, one prop- 
agating into the slower shell (forward shock) and another 
one into the faster shell (reverse shock). The interaction 
of the blobs is modeled as a collision of two homogeneous 
shells flSpada et aL^ OOl JTanihata e t al. 1 2003). and also 
directly simulated (MAMB04). These models shells are as- 
sumed to have sharp edges (density discontinuities) sepa- 
rating them from the background fluid. 

The present work relies on the physical model devel- 
oped by MAMB04. In that model the comoving magnetic 
field in the region behind the shocks is assumed to be 
randomly oriented in space. Its energy density is propor- 
tional to the internal energy density of the shocked fluid. 
The non-thermal particles simulated within the model are 
accelerated at the shock front to very high Lorentz factors 
whereby they radiate synchrotron radiation. 

The observed light curve from a two-shell collision de- 
pends on the position along the jet at which the interac- 
tion occurs as well as on the details of the evolution. As we 
discussed in MAMB04, the hydrodynamic evolution of the 
shells begins already before the collision time predicted by 
the analytic models (e.g. |Spada et al. | 2001) and its most 
important features are (MAMB04): 

— The front (with respect to the direction of motion) 
discontinuity of each shell decays into a bow shock, a 
contact discontinuity and a reverse shock. 



— The back discontinuity of each shell develops into two 
rarefactions. The first one connects the still unper- 
turbed material inside the shell with a contact discon- 
tinuity separating shell matter and external medium. 
The second rarefaction connects the contact disconti- 
nuity with the external medium. 

The pre-collision evolution is quantitatively similar in- 
dependent of whether the edges of the colliding shells are 
modeled as sharp discontinuities or by a smooth transition 
layer. The pre-collision evolution of the shells is important 
as the internal energy of the front part of the shells rises 
due to the supersonic propagation of the shells. Hence, 
when they actually start colliding the conditions inside 
the shells are different from those at ejection from the cen- 
tral engine, i.e., the conditions at the contact surfaces of 
the two shells are different from those of the initial shells. 
Once the shells start to collide the evolution consists of 
three basic stages: 

— The pressure in the region between the two shells rises 
rapidly (compared to the light crossing time of the 
shell). The front and the reverse shock propagate into 
the slower and faster shell, respectively. 

— The shocks heat up the cold shell material. While the 
heated region expands the pressure decreases slightly. 

— The front shock breaks out of the slower shell, enters 
the less dense external medium, and thus accelerates. 
The pressure decreases more rapidly during this stage 
because of the faster expansion. 

3. Numerical simulations 

In our previous paper we simulated shell collisions in two 
spatial dimensions. The results showed that the lateral ex- 
pansion is negligible during the collision of aligned shells 
(MAMB04, Fig. 4). Thus, all essential features of aligned 
shell interactions can be captured using ID simulations. 
This allowed us to compute seven high-resolution ID mod- 
els all having the same initial shell velocities and the same 
shell geometries. The models include one where the rest 
mass densities of both shells are identical (S10-F10), two 
where both shells have approximately the same conserved 
mass (SI0-F07, SI4-F10), and one where both shells pos- 
sess the same kinetic energy (SI0-F05). The latter three 
models were also simulated by exchanging the initial rest 
mass density of the colliding shells (Tab.^). 

The numerical grid consists of 10 4 zones covering a 
physical domain with a length of 5 I0 15 cm. A re-mapping 
technique (see MAMB04 for details) was applied in or- 
der to be able to follow the evolution of the two shells 
initially separated by I0 14 cm up to distances of I0 17 cm 
from the AGN engine. The re-mapping technique allowed 
us to resolve spatial scales from 5 10 11 cm to 10 cm. 

The energy distribution of the non-thermal electrons 
is covered by 64 energy bins, and the number of frequen- 
cies at which we compute the synchrotron radiation is 25, 
logarithmically spanning the frequency range from 10 16 
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Table 1. Overview of the shell properties (densities p\ and p2, Lorentz factors Ti and T2, and inertial masses M.\ 
and M2) of the seven ID models. The initial thickness of the shells is 10 14 cm, their volume is V, and their initial 
separation is Dq = 10 14 cm. The density of the external medium p cx t = 10~ 23 g cm~ 3 , and its Lorentz factor F ox t = 2.9. 
Initially p/pc 2 = 10~ 4 everywhere. The parameters of the type-E shock acceleration model (see MAMB04 for details) 
are: a c = 10 -2 , 7 m i n = 30, r\ — l 10 3 . The last two columns give the Lorentz factor of the reverse (r^ lm ) and forward 
(rp m ) shocks, computed from the hydrodynamic simulations. 



to 10 19 Hz. We use the type-E shock acceleration model of 
MAMB04. 

The integration of the conservation laws of relativistic 
hydrodynamics and the computation of the synchrotron 
radiation are performed using RGENESIS, which is an 
extension of the code GENESIS ( |Aloy et al. 1 1999). The 
handling of the 11011 thermal particles within RGENESIS 
is described in detail in MAMB04. 

3.1. Hydrodynamic evolution 

The main properties of the hydrodynamic evolution of our 
models are summarized in Figs. ^ and |21 showing the initial 
models (top row), the formation of the internal shocks 
(middle row), and the final state (bottom row). 

The evolution depends on the density of the shells in 
the following way: the higher the density of the faster shell, 
the higher is the resulting density peak (middle rows of 
Figs. n and 0), and the larger are the density variations 
within the merged shell (bottom rows in Figs.^ an d 0). 
The shock Lorentz factors, as computed from the simula- 
tions (r| lm , Tj 1111 ), increase with the rest mass density of 
the faster shell. 

Kino et al. (2004) obtain a multi-peaked final den- 
sity distribution, which they attribute to the rarefaction 
waves forming when the internal shocks cross the edges of 
the shells. In our models the rarefaction waves cause the 
Lorentz factor to decrease towards the rear part of the 
merged shell (bottom rows in Figs.nand[2Jl. However, the 
direct comparison with IKino et al. 1 (2004) is difficult, as 
their results are computed in the comoving frame of the 
contact discontinuity between the two shells, whereas our 
results are shown in the laboratory frame attached to the 
central engine. Despite of this difficulty, there is a robust 
result common to every simulation modeling two-shell in- 
teractions, namely the formation of a very inhomogeneous 
merged structure as a result of the non-linear evolution 
of the system. This heterogeneity flaws simple one-zone 
models assuming homogeneous conditions in the emitting 
regions. 



3.2. Energy conversion efficiency 

Assuming an inelastic collision of two shells, the fol- 
lowing conversion efficiency e of bulk kinetic energy 
into internal energy can be derived from the conserva- 
tion of energy and momentum (|Kobayas hi et al.~| 1997; 
|Daigne fc Mochkovitch| 1998): 



1 



(M 1 + M 2 )r m 
M 1 r 1 +M 2 r 2 



(i) 



where Aii — piTiV is the total inertial mass of the shell 
i (i = 1,2), and where the Lorentz factor of the merged 
shell is given by 



r m — \/riF 



MiT 1 +M2T2 



(2) 



■M1T2+M2T1 

The (time independent) efficiency given in Eq. Q cor- 
responds to a situation where two incompressible shells 
collide inelastically and instantaneously. However, as we 
are modeling the shells as compressible fluids, we have to 
consider the temporal evolution of the energy conversion 
efficiency, and hence use the following generalized defini- 
tion 



e(T) 



1 



jT 2 (T,z)p(T,z)dz 
/F2(0,z)p(0,z)dz ' 



(3) 



where p(T, z) and T(T, z) are the rest mass density and 
the Lorentz factor of the fluid at a time T at position z, 
respectively. We point out that the instantaneous values 
of e(T) do not necessarily coincide with those given by 
Eq. because during the hydrodynamic evolution the 
fluid suffers multiple compressions and expansions chang- 
ing the ratio of internal energy and kinetic energy (which 
is assumed to be constant in the simple model). 

Fig-El shows the temporal evolution of conversion effi- 
ciency e(T) for all of our simulated models. Initially, the 
efficiency rises until about 150 ks, and then remains nearly 
constant until 350 ks. The analysis of the hydrodynamic 
evolution shows that this is due to the fact that the faster 
shell is heated due to its interaction with the medium 
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Fig. 1. Hydrodynamic evolution of the shell collision for models S10-F10, S10-F14, S10-F07, and S10-F05. The solid 
lines give the rest mass density in units of p e xt- The dashed lines show the Lorentz factor of the fluid which is moving 
to the right. The origin of the z-axis corresponds to a distance of 0cm (top row), 2.65 10 16 cm (middle row), and 
6.61 10 16 cm (bottom row), respectively. Times given in the leftmost panels of each row are measured in the laboratory 
frame. 




T[ks] 

Fig. 3. Temporal evolution of the conversion efficiency 
defined by Eq. @ for all models. T is the time measured 
in the source frame. 



ahead of it. However, at 150 ks the shell encounters a rar- 
efaction created by the slower shell reducing the amount 



of heating. At 330 ks the internal shocks form and the ef- 
ficiency starts rising again until it reaches a maximum at 
about 1000 ks. The subsequent decrease of the efficiency 
is caused by the expansion of the fluid as the shocks break 
out of and rarefactions form within the merged shell. The 
efficiency is generally the higher the smaller is the rest- 
mass density of the faster shell (Tab.[2J . We note that the 
final values approach the efficiency predicted by the con- 
servation laws (Eq. (JJJ). 

The time evolution of the conversion efficiency 
in our models is qualitatively similar to that of 
Kin o et al. f 2004'). However, due to the fact that these 
authors computed the evolution in the rest frame of the 
(relativistically moving) contact discontinuity, a quantita- 
tive comparison with our results is presently impossible. 

3.3. Synthetic light curves 

The peak intensity of the light curves in both the soft 
(Fig. ^] upper panel; 0.1—1 keV) and the hard energy band 
(Fig. 21 lower panel; 2 — lOkeV) is correlated to the to- 
tal rest- mass of the shells (see Tab. Pi . This correlation is 
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Fig. 2. Same as Fig.HJ but for models S14-F10, S07-F10 and S05-F10. Model S10-F10 is shown once again here to 
allow for an easier comparison. 
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Table 2. Maximum conversion efficiency e max , final con- 
version efficiency £fi n , and theoretically predicted conver- 
sion efficiency 6p rc d (see Eq. QJ) for the seven simulated 
models. 



slightly better in the wider hard energy band. We further 
find that the peak of the flare occurs in the hard energy 
band several hundreds of seconds before the peak in the 
soft band. 

The shape of the flares is qualitatively similar for all 
models indicating that the variation of the density ratio 
between shells mostly tends to change the normalization 
of the flares. Furthermore, the large variations of the den- 
sity (Figs.Q] and [21 bottom rows) and of the specific in- 
ternal energy are not correspondingly do not change the 



shape of the light curves which are rather smooth. This 
implies that from the shape of the flare alone, it is very 
hard to infer any information about the rest-mass den- 
sity or the specific internal energy of the emitting regions. 
This restriction should be taken into account by any ana- 
lytic model aiming to infer physical parameters from light 
curves. 



3.4. Sensitivity to the initial conditions 

The light curves generated from our models are rather 
insensitive to both the choice of density and pressure of 
the external medium provided they are sufficiently small 
to be inertially negligible. If the density of the underlying 
medium is comparable to that of the shells the evolution 
will be quantitatively different (but qualitatively, the same 
Riemann structure will develop from every edge of every 
shell). We have limited the Lorentz factor of the external 
medium to be smaller than that of the slowest shell. Its 
value changes only the depth of the rarefaction trailing 
both shells (the smaller the velocity of the jet, the smaller 
the pressure and the density reached behind every shell). 
However, this will not change the light curves of our mod- 
els as from rarefactions there is no emission (because only 
shocked regions emit in our model). The complementary 
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Fig. 4. Soft (upper panel; observed photon energies be- 
tween 0.1 and IkeV) and hard (lower panel; 2-10 keV) 
light curves for all of our seven models. A correlation be- 
tween the total rest- mass of the shells and the peak photon 
counts is seen. 



situation, i.e., one or both shells being slower than the jet, 
we have not considered. This situation will correspond to 
a couple of shells being slower than the underlying, light 
jet. In such a case, the Riemann structure will be qualita- 
tively different. The rear edges of the shells will be heated 
because of the development of shocks while the pull of the 
external medium will create rarefactions from the forward 
edges of every shell. Nevertheless, this interesting situa- 
tion falls beyond the scope of the current work and will 
be considered elsewhere. We point out that it is reasonable 
to expect a low-density external medium, since the shells 
will clean up the medium through which they propagate 
(note that the state left behind every shell is less dense 
than that set up initially) . This rarefication process of the 
external medium is a robust feature associated with the 
shell motion (see MAMB04). 



4. Analytic modeling of flares 

In this section we discuss some analytic approaches to 
the modeling of flares resulting from internal shocks, and 
most importantly which physical parameters can be ob- 
tained by analyzing the shape and duration of flares. We 
begin with general considerations of the observations of 
distant unresolved sources, and then provide an analytic 
model particularly suited to study the flares resulting from 
collisions of shells in a relativistic jet. 

4.1. Computation of the radiation observed by a 
distant observer 

As the simplest case, we assume that we have a source of 
radiation moving along the line-of-sight (LOS) towards 
a distant observer. We introduce two inertial frames: the 
source frame and the observer frame. The source frame is 
attached to the central energy source. In this frame the 
physics of the source is best described because time and 
space are independent coordinates. In the source frame, 
we denote the position along the LOS by z, and the time 
by T. The observer frame is located at a distance z — z Q b s 
along the LOS. Except for a cosmological redshift factor, 
this frame is equivalent to one attached to the Earth. In 
this frame we have chosen two different coordinate sys- 
tems, the (z bs, £obs)-system and the (x, y)-system. The 
second coordinate system is best suited for describing the 
structure of the observations in space time because the 
coordinates associated with this frame mix the spatial (z) 
and temporal (T) coordinates of the source frame in the 
following way: 



x = 



y 



cT-z 
cT + z 

~7T 



(4) 



(5) 



where c is the speed of light. 

A process which takes place at time T and position z 
in the source frame is observed in the observer frame at a 
time 



^obs — T + 



^■obs Z 



\Plx ^ ^obs 
C C 



(6) 



i.e., all points with the same x coordinate will be observed 
simultaneously. Hence, the total radiation observed by a 
distant observer at time t is equal to the integral over y 
of the emitted radiation of all points in the source frame 
which have the same x — (ct— z bs)/\/2 coordinate in the 
observer frame (xy-coordinate system), i.e., 



!{t) = / dyj(x,y) 



(7) 



where J(i) is the intensity of the observed radiation at 
time t, and j(x, y) is the emissivity in the observer frame. 
This emissivity is related to the emissivity at the corre- 
sponding time and place in the source frame, i.e., j(x, y) = 
j(z(x,y),T(x,y)). 
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If a point source is moving with a velocity [3c then 
z(T) = zq + (3cT, and its trajectory in ^-coordinates will 
be 



y{x) = - — -x 



V2 



-,zo ■ 



(8) 



1-/3 1-/3 

Hence, an object moving with the speed of light is ob- 
served to propagate parallel to the y axis, while an object 
at rest in the source frame is observed to move parallel to 
the line y = x. 

Now we assume that an observer is located at z b s = 0, 
and that the emissivity of a point source is a function of 
the time T , 



3{z,T) 



I S{z -zq- j3cT) if < T < T c 



end 



(9) 







otherwise 



where Iq, zq, (3 and T on( j are the normalization of the emis- 
sivity, the initial source position, its velocity in units of the 
speed of light, and the duration of its emission, respec- 
tively. In the observer frame (rry-coordinates) the emis- 
sivity is then, using |@J, (0, and the properties of the 
(5-function, 



J v2 

j{x,y) = - — -o 



x(l + (3) + zqV2 
1-/3" T l-f3 

S(y;-x, V2cT cn d-x) 



(10) 




1 if a < x < b 
otherwise 



(11) 



Using now equation (|7|) we can compute the light curve 
produced by a single emitting point source observed by a 
distant observer: 



Ct -Zq (1 - (3)cT cnd 



(12) 



0~ VV2' V2' V2 

It follows from Eq. I|12|l that for the observer the duration 
of the emission is given by 



ttotal = (1 — P)T e nd , 



(13) 



or, in the ultra-relativistic limit by t to ta\ ~ T cnd / (2T 2 ), 
where F = 1/ y/l — /3 2 is the Lorentz factor of the emitting 
point source. 

This shows that the total time (measured in the source 
frame) during which the source emits and its Lorentz fac- 
tor are degenerate for a distant observer, i.e., it is not pos- 
sible to determine both the Lorentz factor and the source 
frame time of an emission process from the duration of the 
observed emission alone. This is a major obstacle for any 
model which attempts to reconstruct, from an observed 
flare, the sizes and the velocities of colliding shells which 
produce that flare. Nevertheless, in the following sections 
we develop an analytic model which can recover some of 
the properties of the colliding shells, taking into account 
this limitation. 



4.2. Geometry of the emitting region 

Motivated by the results of numerical simulations 
(MAMB04) the emitting flare can be represented as a 
horn-shaped region in the observer frame when using the 
^-coordinates (Fig.[5l). The two curves yb(x) and yi{x) 
delimiting the horn-shaped region correspond to the front 
and back edges of the merged shell, respectively. They 
are approximated by second order polynomials passing 
through the origin (Fig.[5J|. This leaves two free param- 
eters for each curve, which are determined by requiring 
that the velocities of the front and back edges be /3g and 
/?bi at T = 0, and /3{ e and /3be at T — T e , respectively. 
T e is the time (measured in the source frame) at which 
the emissivity from the shell interaction drops below the 
observable level. By computing the equations of motion 
(in the observer frame) Zf(T) and Z\,{T) corresponding to 
yt(x) and j/b(x), respectively, one obtains the following 
equations for the front edge, 

n X 1 

K f T 2 r 



yt(x) 



z f {T) =cT{l 



4F?x 



V2T C 



4Fj? 



1-77^2 - (14) 



r 2 
1 r 



1 + K r 



T 

c t1 



and the back edge 




where Ff, Tb, K{, and arc defined as 
F f ^[2(l-/3 fi )]- 1/2 , 
Fb^[2(l-/? bi )r 1/2 , 



1-/3b 

I -Pie 
l-/?bi 



1 



1 . 



(15) 

(16) 
(17) 



(18) 
(19) 

(20) 
(21) 



We point out that only in the ultra-relativistic case 
(Lorentz factors XL 2, j3 close to 1) Tf and Tb approach 
the Lorentz factors corresponding to the initial velocities 
of the front and back edges which in turn can be identi- 
fied with the Lorentz factors of the forward and reverse 
shocks, respectively. However, as we are usually dealing 
with ultra-relativistic velocities in our model, there arises 
no significant error when substituting l|18|l and i|19fl by 
the corresponding Lorentz factors in equations i|14|) - (|17|) . 
In this case Kf and K\, are given by 

4 



K t = ^ - 1 



K h 



r b 



i 



(22) 



(23) 



where Tf e and Tbe are the Lorentz factors of the front and 
back edge at T = T e , respectively. 
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Fig. 5. Light curve (upper panel) and the observer emis- 
sivity distribution in xy-coordinates (lower panel) for a 
horn-shaped segment with parameters I\ = 8, Kb = 4, 
r f = 12, Ki = 3, Ti = 30 ks, T 2 = 100 ks, T c = 200 ks, 
r] = 0.1, and 6 = 1. The curves delimiting the horn-shaped 
region corresponding to the front and back edges of the 
merged shell are denoted by yt(x) and yb(x), respectively. 
The edges of the segment corresponding to T% and T 2 are 
indicated by y\{x) and ys(x), respectively. The vertical 
lines correspond to the points x\t (solid), x\b (dotted), x%t 
(dashed), and x 2 b (dot-dashed), respectively. The num- 
bers inside the horn-shaped area denote regions where 
the emissivity is integrated using different upper and lower 
limits. 

4.3. Light curve of a horn-shaped segment 

In order to compute the light curve emitted by the horn- 
shaped region (Fig.|SJ, we note that significant emission 
occurs only after the actual interaction of the shells starts 
and that this emission only lasts until the internal shocks 
reach the edge of the merged shell which terminates the 
shell interaction. Hence, one needs to calculate only the 
light curve produced by a segment of the horn-shaped 
region denned by two times, T\ (start of the interaction) 
and T 2 (termination of the interaction). We assume an 
emissivity distribution in the source frame of the form 



j(z,T) = j Q 



1 



± 2 ± l 

SpiT^Sfozb^ZiiT)] , 



(24) 



where 9 is a power law index describing the temporal 
change of the emissivity, and r\ is the ratio of the final 
j(T 2 ) and the initial jo = j(Ti) emissivities. The interval 
function S(T; T\, T 2 ) (JUJ appearing in the above expres- 
sion will be dropped in the following discussion. 



Substituting Q and (JSJ into J23J, the emissivity in the 
observer frame (xj/-coordinates) becomes 



j{x,y) 



Jo 



Tf - Tf 



jx + y) 6 



fa - 1) 



(25) 



The segment of the horn-shaped region is bound by the 
two parabolae yt(x) and yb{x) in x-direction (see § 14. 2|) . 
and by the following two straight lines in y-direction 
(Fig.0, 



yi (x) = V2cTi - x , 
y 2 (x) = v2cT 2 - x . 



(26) 
(27) 



The light curve I(x) at a given x-coordinate is ob- 
tained by integrating the emissivity distribution at that 
x-coordinate in y-direction. According to our model the 
integration limits depend on the value of the x-coordinate, 
and three regions can be distinguished (Fig.|5J. In region 1 
(xif < x < xib) the integration extends from yi(x) to 
yi(x), in region 2 (xi b < x < x 2 f) from y h (x) to yt(x), 
and in region 3 (x 2 f < x < x 2 b) from yb(x) to y 2 (x) (see 
Fig.EJ. The coordinates xib, Xif, x 2 b, and x 2 f are given 

by 



Xib 



Xif 



X2b 



X2i 



V2cT x 




1 


in 1+ 




f KbTjTe 


V2cTi 




1 


21* 1 + 


vi 


+ K f T 1 /T e 


V2cT 2 




1 


2rg 1 + 


n/i- 


f K h T 2 /T e 


V2cT 2 




1 


2rf i 4 


Vi 


+ K f T 2 /T e 


we assume 


Xib 


< X21 (see 



(28) 
(29) 
(30) 
(31) 



sion). Thus, xif < x±b < x 2 i < X2b- Defining two auxiliary 
functions 



5 (£,if,r)-4r 2 (r 2 !^ 2 + £), 

and 

1 



-6 _ ^8 



(T§-r)T?)(B-A) + 



(r) - l)(B e+1 - A 6 



1 



(32) 



(33) 



allows us to write the light curve of a segment of the horn- 
shaped region in a very compact way, 

I(x) =j Q V2cT e x 



H 
H 
H 



Ti Ti T 2 

-*-e -*- e -*■ e 



9b(x),g{(x), p-, p,6»,?7 



, ,t 2 t 1 t 2 

9b[X), — , — , jT,0,77 



if xif < x < x 



lb 



(34) 



if xib < x < x 2 f 
if x 2 f < x < x 2b 
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where gt(x) and gb{x) are given by 

The upper panel of Fig. [S] shows an example of a light 
curve computed from Eq. l|3"4*|) for Tb = 8, K\, = 4, Tf = 

12, K { = 3, Tx = 30 ks, T 2 = 100 ks, T c = 200 ks, 77 = 0-1, 
and 6 = 1. 

4.4. Two-phase flare model 

As discussed in §0 and by MAMB04, the interaction of 
two shells consists of three basic phases: a fast rise of the 
pressure, a slow pressure decrease while the internal shocks 
propagate through the shells, and finally the acceleration 
of the front shock as it breaks out of the slower shell. In 
our analytic model we simplify the situation by assum- 
ing only two phases, and hence call it a two-phase flare 
model. For the first phase (0 < T < TQ we assume that 
the initial and final values of the emissivity are jo and 
respectively. In the second phase (7] < T < T c ) the 
emissivity is assumed to abruptly jump to a value 772 jo , 
and then to decrease to zero. The abrupt jump is moti- 
vated by the fact that when the forward shock, created 
by the interaction of the shells, is about to break out of 
the slower shell, it encounters the reverse shock produced 
by the pre-collision hydrodynamic evolution of its leading 
edge. The interaction of these two shocks causes a fast rise 
both of the pressure and the density, i.e., 772 > r)x- 

The light curve I(x) of the two-phase flare model is 
then computed as follows: we assume that the emissivity 
j(T) has the form 




x/c [ks] 

Fig. 6. Two-phase flare model with Tb = 5.5, K\, = 1.1, 
r f = 6, Kf = 4, Ti = 50 ks, T c = 100 ks, 6 X = 1, 6 2 = 1 
and 771 = 100, and 772 = 300. The upper panel shows the 
normalized light curve, and the lower one the emissivity. 
The vertical lines denote the coordinates X\ (solid), X2 
(dotted), x 3 (dashed), and x± (dot-dashed line) which di- 
vide the emissivity distribution into four distinct regions 
labeled by numbers 1 to 4. The white thick horizontal line 
separates the two phases of the temporal evolution of the 
emissivity. 



(?7l-l) if0<T<Ti 



3{T)=h{ 



ft>2 X 



T" 



(35) 



if Ti <T <T C 



and that the light curve in both phases (i.e., in both sub- 
segments of the horn-shaped emission region; see (Fig. 01 
has the form of the light curve I(x) 134|) with common pa- 
rameters Tb, -Kb, Tf, and K{, but with the sub-segment 
specific parameters according to the following substitu- 
tions: 



segment 1: 
segment 2: 
77 -> . 



jo, Tx -> 0, T 2 
Jo -> mh, Tx 



Ti, 6- 
T h T 2 



>T C ,6^ 



In order to obtain the light curve T{x) of the two-phase 
flare model (again by integrating the emissivity distribu- 
tion at each x-coordinatc in y-dircction) the horn shaped 



emissivity region is divided into four distinct regions de- 
fined by (Fig.EJ) 



xx 



X-2 



3'3 



V2cTi 



2T 2 { 
V2cTi 



2TI 

V2cT e 



2T 2 { 
V2cT e 



2TI 



l + y/l+ KfTi/Tc 



1 + ^l + K h Ti/T c 



1 



1 



K f ' 



1 + VTTkI 



(36) 
(37) 
(38) 
(39) 
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Fig. 7. Same as Fig.® but for 772 = 101. 

Then, using the auxiliary functions H , g b and <?f defined 
in § 14.31 one can write the light curve of the two-phase 
model in the form 

l{x) = joV2cT e x 



H 



H 



g b (x),g f (x),0, 7 jr,0 1 ,T) 1 



if < x < 



V2H 



V2H 



Ti Ti 



g h (x),g { (x), j^-, 1,6*2,0 



9b{x), 1, j 1, 6*2,0 



if Xi < X < X2 



(40) 



if x-i < x < X3 



if X3 < x < X4 



Figs.|ni and [7] show two examples of analytic two-phase 
light curves, which only differ in the value of 772 with 7/2 = 
300 in Fig. 0and 772 = 101 in Fig.Q respectively. 

4.5. Normalized two- phase model 

As pointed out in S I4.1I all observed time scales are func- 
tions of T/r 2 , where T is a typical time scale in the source 
frame, and T is a typical Lorentz factor at which the emit- 
ting region is moving. In order to make our model "aware" 
of this, we introduce the normalized coordinate 



X4 



op2 

b -(1 + 71 + ^): 



V2cT c 



(41) 



and define 

7b 



r f 



and 



(42) 



(43) 



Thereby we can specify a normalized (both in intensity 
and in time) flare model 

' H[h h (O,h { (O,0,tu6uVi] if < C < Ci 

fl> b (O^i,0,iiA,7 7l ] + 
F{0=Fo{ V2 H[ti,h f (O,ti,l,e 2 ,0] lfa<C - 6 (44) 
V2H[h h (0,h ( (0, U, 1, 6-2, 0] if 6 < £ < 6 



[ 772fr[/i b (0, Mi, i,02,o] 



if & < £ < 1 



where the constant Tq is chosen such that the maximum 
value of •/-"(£) is 1. The points £j, i = 1, 2, 3 are given by 



6 



2+ i + vm?b" 



1 + VT+lMi 



1 + VT+Kb 



T + VI + K h U 



2 i + yi + A b 

7b 1 + ym^ 



(45) 
(46) 
(47) 



and the auxiliary functions /ib(£) and ft.f(£) are obtained 
from <7 b (x) and gi(x), respectively: 



MO 
and 

MO : 



Ab£ 2 



+ 



2£ 



(1 + vi + A-,,) 2 1 + vm?b" 



A f £ 2 



2i 



7^(1 + VT+Ab) 2 7b (! + vT+Ab) 



(48) 



(49) 



The normalized two-phase model has the advantage 
that the absolute intensity is no longer a parameter, and 
that it has embedded the intrinsic degeneracies of T/T 2 
which must be present in the model due to the relativis- 
tic time contraction. Indeed, the normalized two-phase 
model uses two parameter less in order to fit the light 
curves. Therefore, we find that the fitting procedure is 
always convergent. 

5. Fits of the analytic model to synthetic light 
curves 

The final goal of our analytic modeling of the light curve 
of a flare is to extract unobservable physical parameters of 
the emitting region. As a first step, we have validated our 
analytic model by comparing the parameters extracted 
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from fits 1 of synthetic light curves (computed from our 
hydrodynamic models) to the parameters that we can di- 
rectly obtain from the same hydrodynamic models. The 
second step will be to apply our validated analytic model 
to real observations, which is beyond the scope of this 
paper. 

Table[3]shows the results of the fit of the analytic model 
to the soft (0.1-1 keV) light curves of all seven models. 
First, we observe that the ratios 7 b :lm = r b lm /Tf lm ) and 
7b = Tb/Tf are very similar, but not exactly the same, 
i.e., there exists a correlation between the values obtained 
from the simulations and the fit. Second, the small values 
of the parameter K\, indicate a lack of acceleration of the 
reverse shock, i.e., the reverse shock moves at almost con- 
stant velocity until it breaks out of the faster shell. This is 
clearly seen in Fig-El where the trajectory (world-line) of 
the reverse shock (full thick line) is a straight line between 
points a and d, which correspond to the start of the in- 
teraction and the time and place where the reverse shock 
breaks out of the faster shell, respectively. Third, the val- 
ues of the parameter K{ imply that the forward shock 
accelerates, though its Lorentz factor does not increase by 
more than about 30% (model S10-F14). From the simu- 
lations, we see that the forward shock (dashed thick line 
in Fig.|nj) moves at approximately constant velocity be- 
tween points a and b the latter point corresponding to 
the time and place where the forward shock catches up 
with the reverse shock emerging from the leading edge of 
the slower shell - see Sec. l4.4fl . Afterwards, beyond point 
b, the forward shock speeds up, because it breaks out of 
the slower shell. Therefore, in the time interval between 
the formation of the forward shock and the time when the 
emissivity associated to it ceases to contribute to the light 
curve, there is a net acceleration of this shock. This is also 
properly captured by our analytic model. Fourth, on the 
one hand, considering the models having the same inertial 
mass of the slower shell (rows 1-4 in Tab. 03 the time 
necessary for the forward shock to break out of that shell, 
relative to the total flare duration in the laboratory frame, 
decreases with increasing inertial mass of the faster shell. 
On the other hand, we find a similar correlation consider- 
ing those models where the inertial mass of the faster shell 
is fixed (rows 1, 5, 6 and 7), i.e., ib is decreasing as the 
inertial mass of the slower shell increases (see Sec. l7.1|) . 

6. Space time emissivity distribution 

From Fig.El which shows the soft light curve (upper panel) 
and a grey contour plot of the emissivity distribution in 
the observer frame using the xy-coordinates (lower panel) 
for model S05-F10, one sees that the emissivity distri- 
bution has the shape expected from the analytic model, 
i.e., a horn-like shape starting at point a where the in- 
ternal shocks form. The region then widens as the shocks 

1 Due to the complex nature of the multidimensional pa- 
rameter space of our analytic model (Eq. I44H - ). we have used 
genetic algorithms to fit it to the normalized light curves of 
our simulations. 




Fig. 8. Soft light curve of model S05-F10 (dotted line), 
and the best fit of the normalized analytic model (full 
line). 

propagate through the shells. The bulk of the radiation is 
produced relatively early in the evolution (up to point c 
in Fig.El which corresponds to the moment of maximum 
emissivity), i.e., the subsequent formation of rarefaction 
waves does not have any influence on the observed light 
curve (see, however, S 17.11 for a discussion of the possi- 
ble influence on later collisions). We point out that there 
might be a difference between the time (in the observer 
frame) when the light curve has a maximum (point e), 
and the point when the emissivity has a maximum (point 
c) . This difference arises because the total observed inten- 
sity at a given time is the result of an integration of the 
emissivity over the emitting region. Thus, the maximum 
observed intensity can result from a region where the emis- 
sivity is not maximal, but whose emitting volume is suffi- 
ciently large. In case of model S05-F10 (Fig.EJ the maxi- 
mal observed intensity is delayed with respect to the max- 
imum of the emissivity. This influences the shape of the 
light curve, because prior to the maximum, at i bs ~ 6ks, 
we see a kink in the light curve where its slope steepens. 
When the emissivity maximum (in the observer frame) 
happens close to, or even after the intensity maximum, no 
such kink will appear (see Figs .0 and [7] which illustrate 
precisely these two cases). 

7. Discussion 

We have performed a detailed study of two-shell colli- 
sions in the framework of the internal shock model using 
relativistic hydrodynamic simulations. The non-thermal 
radiation emitted by shocked plasma is consistently cou- 
pled to the hydrodynamics (including the radiation back- 
reaction on the dynamics) in our numerical scheme. Our 
findings can be applied to extract physical parameters of 
the emitting regions of blazar flares. However, the method 
of extracting physical parameters out of synchrotron light 
curves might as well be applied to other astrophysical 
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Table 3. Results of the fit of the normalized analytic model to the soft (0.1-1 keV) light curves of all simulated 
models. 7§ lm is the ratio r|, lm /rf lm computed from the models (see Tab.^). The rest of the columns provide the other 
parameters obtained by fitting the normalized analytic model. 
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Fig. 9. The distribution of the emissivity in the observer frame in xy-coordinates in the soft photon band (lower 
panel), and the corresponding soft light curve (upper panel) for model S05-F10. The grey contours show the emissivity 
distribution on a linear scale. The full, dotted, dashed, and dot-dashed lines denote positions of the back edge of the 
slower shell, the front edge of the slower shell, the back edge of the faster shell, and the front edge of the faster shell, 
respectively. The full thick and dashed thick lines show the trajectories of the reverse and forward shock, respectively. 



scenarios where internal shocks occur (e.g., gamma-ray 
bursts). 

7.1. Hydrodynamic evolution of the shell interaction 

The properties of the interaction common to all our mod- 
els are the existence of three main phases of evolution 



(Figs.Q] and [2J : (i) the pre-collision phase where the 
front edges of the shells pile up matter from the exter- 
nal medium and are heated due to their interaction with 
the latter, (ii) the collision phase where the internal shocks 
form and propagate through the shells, and (iii) the post- 
collision or rarefaction phase where the reverse shock has 
broken out of the faster shell creating rarefaction waves 
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which expand the merged shell and cause its structure to 
become multi-peaked. Although this multi-peaked struc- 
ture does not have any influence on the observed flare, one 
has to bear in mind that any shell colliding subsequently 
with the merged shell will encounter a non-uniform struc- 
ture, which will probably cause the forward shock, orig- 
inating after the new shell collision, to undergo a much 
more complicated set of non-linear wave interactions than 
in the case of initially uniform shells. This will yield a more 
complex light curve profile. We also point out that the in- 
sensitivity of the light curve profiles to the large rest-mass 
density and specific energy variations in the merged shell, 
may probably prevent one from extracting from any sim- 
ple one-zone model accurate information about these hy- 
drodynamic variables from fits of the observed flare light 
curves. This lack of sensitivity to the rest-mass density 
profile comes from the fact that most of the steep den- 
sity gradients in the merged shell are either weak shocks 
or no shocks at all. Hence, no significant emission can be 
originated in such steep density variations and the result- 
ing light curve remains uninfluenced. Furthermore, when 
assuming an approximately uniform proportionality be- 
tween the rest-mass energy and the magnetic field energy 
(i.e., the ratio B 2 /p being approximately uniform) accu- 
rate values of the magnetic field strength are hardly to be 
inferred from a flare light curve, too. 

We have also computed the instantaneous efficiency of 
conversion of kinetic energy into internal energy of the 
fluid. As can be seen from Fig. [31 the efficiency rises ini- 
tially, then remains constant for a period of time during 
which the front edge of the faster shell propagates through 
the rarefaction caused by the back edge of the slower shell, 
and then rises again when the internal shocks form. The 
efficiency during the period when most of the radiation is 
observed is larger than the analytically predicted value. 
Afterwards, it decreases and reaches the analytically pre- 
dicted values. 

For a fixed inertial mass of the slower shell, the Lorentz 
factor of the reverse and forward shocks grows as the in- 
ertial mass of the faster shell increases. Conversely, fixing 
the inertial mass of the faster shell, the Lorentz factor of 
the shocks decreases with increasing inertial mass of the 
slower shell. These correlations are expected, because in- 
creasing the inertial mass of the slower shell any incoming 
(faster) shell will supply less inertia (relative to that of 
the slower shell), i.e., the two internal shocks will propa- 
gate at smaller velocities. On the other hand, increasing 
the inertia of the faster shell will result in larger shock 
speeds, because the target shell does not decelerate the 
faster shell so efficiently. 

The parameter t\, which is approximately the ratio of 
the time it takes the forward and reverse shocks to cross 
the slower and faster shell, respectively, depends on the 
inertial masses of the shells, too. We find that U decreases 
with increasing relative inertia of the faster shell (with 
respect to the slower one), because with increased inertia it 
is much easier for the forward shock to propagate through 
the slower shell and exiting from its front edge earlier than 



the reverse shock, which propagates into the heavier, faster 
shell. 

7.2. Properties of the analytic model 

As described in §0J we have introduced the observer-xy 
coordinate system in order to simplify the relationship be- 
tween the observed radiation and the space time evolution 
of the emissivity. By making assumptions about the prop- 
erties of the emissivity distribution, we are able to con- 
struct an analytic model which depends on parameters of 
the internal shocks, namely the shock Lorentz factors and 
the shock propagation time scales. Our model can be ap- 
plied to non-relativistic as well as ultra-relativistic shocks, 
and even allows for a (constant) acceleration of the shocks. 
For the moment, it does not take into account any spectral 
information, but we intend to improve our model in order 
to be able to interpret multi-frequency observations which 
may provide further restrictions on the physical parame- 
ters of the emitting regions. We have validated the analytic 
model by comparing its fitted parameters (using synthetic 
light curves) with their values obtained directly from the 
simulations which provide the synthetic light curves. 

7.3. Synthetic flares and the analytic model fitting 

We find that all flares, once normalized, look qualitatively 
the same. This suggests that the exact shape of the flare 
depends mostly on the shock Lorentz factors and shock 
propagation time scales while it is rather insensitive to 
the distribution of the rest-mass density or specific energy 
within the emitting region (see S I7.1[) . We have fitted the 
analytic model (using a simple genetic algorithm) to syn- 
thetic light curves and have achieved satisfactory results. 
The deviations in the parameters between the simulation 
and the fit (Tab.|2J are due to the simplifications of the 
analytic model, e.g., the assumption that the emissivity is 
constant between internal shocks (this is not exactly the 
case, as can be seen in Fig.EJ). However, the general shape 
of the emissivity distribution (Fig.E|) does match our as- 
sumptions. This implies that one might fit the model pa- 
rameters using the observed X-ray light curves of blazar 
flares (e.g., Mrk421) and recover shock parameters, like 
e.g., the ratio of Lorentz factors between the reverse and 
forward shock, the propagation time scales through the 
shells, and the kinematic evolution (acceleration) of the 
shock fronts. However, here the quality of the observa- 
tions might play a crucial role, since the analytic model 
is very sensitive to changes in the flare shape, as can be 
seen from Figs. El and 

We have produced synthetic observations in soft and 
hard energy bands (Fig.0J. A correlation between the 
peak photon counts and the initial rest mass density of 
the shells has been found. Currently, we are performing 
a more detailed parameter study aiming to confirm this 
suspected correlation for cases of different shell Lorentz 
factors and densities. The peak in the hard band is ob- 
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served several hundreds of seconds before the one in the 
soft band. However, the exact value of the lag between 
the hard and soft peaks will be different in the case of 
varying shock strengths and velocities (which strongly de- 
pend on the shell properties). Additionally, the possible 
acceleration of particles inside the emitting region is not 
included in our model. Instead, particles are injected into 
the plasma with a given energy distribution, as acceler- 
ation timescales are expected to be shorter than the dy- 
namical timescale. Note, however, that acceleration may 
be the origin of some of the features (sign of the de- 
lay) observed in temporal evolution of these sources (e.g., 
|Kirk, Rieger fc Mastichiadis| 1998). 

Observationally, the lags between the peaks in the 
soft and hard bands are still a matter of controversial 
debate. The emission of the soft X-rays can be well 
correlated with that of the hard X-rays lagging it by 
3-4 ks IjTakahashi et al~l 1996, 2000 , |Zhang et al. 1 1999, 
IMalizia et aTl 2000, IKataoka et"ai~l 2000. IFossati et al. I 
2000). However, significant lags of both signs were de- 
tected from several flares ifTanihata et al. 1 2001). The lags 
between the soft and hard peaks may depend on the as- 
sumed value of the parameter olb = (B 2 / 8ir) / (pe) which 
sets the strength of the magnetic field in the emitting re- 
gion. As in MAMB04 this parameter has been chosen such 
that the field strength resulting from the hydrodynamic 
evolution yields values of ~ 0.1 G (the value of as is the 
same for all our models). This field strength lies in the 
ballpark of the values inferred from detailed fits of SED 
curves of several sources (see, e.g.. ITanihat a et al. 1 2001). 
We point out that even for the same source there is no con- 
sensus on the exact value of the magnetic field strength. 
For example. ITanihata et aT~l f2001) report B ~ 0.14 G 
for PKS 2155-304 , while is inferred B ~ 1.5 G from 
|Chiappetti et al. 1 (1999) for the same source (taking the 
Doppler factor to be 5 = 28: ITanihata et al. 1 2001). 

Observations of PKS 2155-304 lEdelson eTaFI (2001) 
using XMM — Newtonsuggest that previous claims of 
soft lags with time scales of ~ hours might be an arti- 
fact of the periodic interruptions of the low-Earth orbits 
of the satellites every ~ 1.6 hours. This claim was ques- 
tioned by |Zhang et al. | (2004) who show that, although 
periodic gaps introduce larger uncertainties than present 
in evenly sampled data, lags on time scales of hours cannot 
be the result of periodic gaps. Large flares with time scales 
of ~ 1 day were detected with temporal lags of less than 
1.5 hours between X-ray and TeV energies (for Mrk421 
see ITakahashi et al. 1 2000) . 

Recently, very detailed studies of Mrk421 with 
XMM - jVew ton l|Brinkmann et all 2003, 2005; 
IRavasio et al. 1 2004) clearly show that the source 
exhibits stronger variability when it is brighter, that the 
cross correlation function appears to change "continu- 
ously" , and depending on the length of the observing 
window and the actual activity state of the source they 
find periods with positive, negative or no lags, but also 
periods of weak correlations between the soft and hard 
energy bands. Here we propose that the apparent sign 



of the lags is the result of the particular arrangement of 
a number of emitting regions, i.e., the resulting spectral 
properties will depend strongly on the time scales of their 
emission, as well as their distances and relative velocities. 
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